home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / STATMISC.ZIP / TETRA.FOR < prev    next >
Text File  |  1985-12-31  |  5KB  |  171 lines

  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE TETRA
  5. C
  6. C     PURPOSE
  7. C        COMPUTE A TETRACHORIC CORRELATION COEFFICIENT BETWEEN TWO
  8. C        VARIABLES WHERE DATA IN BOTH VARIABLES HAVE BEEN REDUCED
  9. C        ARTIFICIALLY TO TWO CATEGORIES.
  10. C
  11. C     USAGE
  12. C        CALL TETRA (N,U,V,HU,HV,R,RS,IE)
  13. C
  14. C     DESCRIPTION OF PARAMETERS
  15. C        N  - NUMBER OF OBSERVATIONS
  16. C        U  - INPUT VECTOR OF LENGTH N CONTAINING THE FIRST VARIABLE
  17. C         REDUCED TO TWO CATEGORIES
  18. C        V  - INPUT VECTOR OF LENGTH N CONTAINING THE SECOND VARIABLE
  19. C         REDUCED TO TWO CATEGORIES
  20. C        HU - INPUT NUMERICAL CODE INDICATING THE HIGHER CATEGORY OF
  21. C         THE FIRST VARIABLE.  IF ANY VALUE OF VARIABLE U IS
  22. C         EQUAL TO OR GREATER THAN HU, IT WILL BE CLASSIFIED AS
  23. C         THE HIGHER CATEGORY, OTHERWISE AS THE LOWER CATEGORY.
  24. C        HV - SAME AS HU EXCEPT THAT HV IS FOR THE SECOND VARIABLE.
  25. C        R  - TETRACHORIC CORRELATION COMPUTED
  26. C        RS - STANDARD ERROR OF TETRACHORIC CORRELATION COMPUTED
  27. C        IE - ERROR CODE
  28. C         0 - NO ERROR
  29. C         1 - UNABLE TO COMPUTE A TETRACHORIC CORRELATION DUE TO
  30. C             THE FACT THAT AT LEAST ONE CELL SHOWS ZERO FRE-
  31. C             QUENCY IN THE 2X2 CONTINGENCY TABLE CONSTRUCTED
  32. C             FROM INPUT DATA.  IN THIS CASE, R AND RS ARE SET
  33. C             TO 10**75.  (SEE GUILFORD, 1956)
  34. C         2 - THE ROOT SOLVER GIVES MULTIPLE ROOTS, OR NO ROOTS,
  35. C             R, IN THE INTERVAL (-1,1) INCLUSIVE. R AND RS ARE
  36. C             SET TO 10**75.
  37. C         3 - UNABLE TO COMPUTE A SATISFACTORY VALUE OF TETRA-
  38. C             CHORIC CORRELATION USING NEWTON-RAPHSON METHOD OF
  39. C             APPROXIMATION TO THE ROOT OF THE EQUATION.  R AND
  40. C             RS ARE SET TO 10**75.  SEE SUBROUTINE POLRT ERROR
  41. C             INDICATORS.
  42. C         4 - HIGH ORDER COEFFICIENT OF THE POLYNOMIAL IS ZERO.
  43. C             SEE SUBROUTINE POLRT ERROR INDICATORS.
  44. C
  45. C     REMARKS
  46. C        VALUES OF VARIABLES U AND V MUST BE NUMERICAL, AND
  47. C        ALPHABETIC AND SPECIAL CHARACTERS MUST NOT BE USED.
  48. C        FOR A DEPENDABLE RESULT FOR TETRACHORIC CORRELATION,
  49. C        IT IS RECOMMENDED THAT N BE AT LEAST 200 OR GREATER.
  50. C
  51. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  52. C        NDTRI
  53. C        POLRT--THIS POLYNOMIAL ROOT ROUTINE WAS SELECTED BECAUSE OF
  54. C           ITS SMALL STORAGE REQUIREMENT.  OTHER SSP ROUTINES
  55. C           WHICH COULD REPLACE POLRT ARE PRQD AND PRBM.  THEIR
  56. C           USE WOULD REQUIRE MODIFICATION OF TETRA.
  57. C
  58. C     METHOD
  59. C        REFER TO J. P. GUILFORD, 'FUNDAMENTAL STATISTICS IN PSYCHO-
  60. C        LOGY AND EDUCATION', MCGRAW-HILL, NEW YORK, 1956, CHAPTER 13
  61. C        AND W. P. ELDERTON, 'FREQUENCY CURVES AND CORRELATION' 4-TH
  62. C        ED., CAMBRIDGE UNIVERSITY PRESS, 1953, CHAPTER 9.
  63. C
  64. C     ..................................................................
  65. C
  66.       SUBROUTINE TETRA (N,U,V,HU,HV,R,RS,IE)
  67. C
  68.       DIMENSION XCOF(8),COF(8),ROOTR(7),ROOTI(7)
  69.       DIMENSION U(1),V(1)
  70.       DOUBLE PRECISION X31,X32,X312,X322
  71. C
  72. C     CONSTRUCT A 2X2 CONTINGENCY TABLE
  73. C
  74.       A=0.0
  75.       B=0.0
  76.       C=0.0
  77.       D=0.0
  78.       DO 40 I=1,N
  79.       IF(U(I)-HU) 10, 25, 25
  80.    10 IF(V(I)-HV) 15, 20, 20
  81.    15 D=D+1.0
  82.       GO TO 40
  83.    20 B=B+1.0
  84.       GO TO 40
  85.    25 IF(V(I)-HV) 30, 35, 35
  86.    30 C=C+1.0
  87.       GO TO 40
  88.    35 A=A+1.0
  89.    40 CONTINUE
  90. C
  91. C     TEST WHETHER ANY CELL IN THE CONTINGENCY TABLE IS ZERO.
  92. C     IF SO, RETURN TO THE CALLING ROUTINE WITH R=0.0 AND IE=1.
  93. C
  94.       IE=0
  95.       IF(A) 60, 60, 45
  96.    45 IF(B) 60, 60, 50
  97.    50 IF(C) 60, 60, 55
  98.    55 IF(D) 60, 60, 70
  99.    60 IE=1
  100.       GO TO 86
  101. C
  102. C     COMPUTE P1, Q1, P2, AND Q2
  103. C
  104.    70 FN=N
  105.       P1=(A+C)/FN
  106.       Q1=(B+D)/FN
  107.       P2=(A+B)/FN
  108.       Q2=(C+D)/FN
  109. C
  110. C     FIND THE STANDARD NORMAL DEVIATES AT Q1 AND Q2, AND THE
  111. C     ORDINATES AT THOSE POINTS
  112. C
  113.       CALL NDTRI (Q1,X1,Y1,ER)
  114.       CALL NDTRI (Q2,X2,Y2,ER)
  115. C
  116. C     COMPUTE THE TETRACHORIC CORRELATION COEFFICIENT
  117. C
  118.       IF(X1) 76, 72, 76
  119.    72 IF(X2) 76, 74, 76
  120.    74 R=0.0
  121.       GO TO 90
  122.    76 XCOF(1)=-((A*D-B*C)/(Y1*Y2*FN*FN))
  123.       XCOF(2)=1.0
  124.       XCOF(3)=X1*X2/2.0
  125.       XCOF(4)=(X1*X1-1.0)*(X2*X2-1.0)/6.0
  126.       X31=DBLE(X1)
  127.       X32=DBLE(X2)
  128.       X312=X31**2
  129.       X322=X32**2
  130.       XCOF(5)=SNGL(X31*(X312-3.0D0)*X32*(X322-3.0D0)/24.0D0)
  131.       XCOF(6)=SNGL((X312*(X312-6.0D0)+3.0D0)*(X322*(X322-6.0D0)+3.0D0)
  132.      1          /120.0D0)
  133.       XCOF(7)=SNGL(X31*(X312*(X312-10.0D0)+15.0D0)*X32*(X322*(X322-10.0
  134.      1          D0)+15.0D0)/720.0D0)
  135.       XCOF(8)=SNGL((((X312-15.0D0)*X312+45.0D0)*X312-15.0D0)*(((X322-
  136.      1          15.0D0)*X322+45.0D0)*X322-15.0D0)/5040.0D0)
  137. C
  138.       CALL POLRT (XCOF,COF,7,ROOTR,ROOTI,IER)
  139. C
  140.       J=0
  141.       IF(IER) 78, 78, 84
  142.    78 DO 82 I=1,7
  143.       IF(ABS(ROOTI(I))-.5*ABS(ROOTR(I))*1.0E-6)79,79,82
  144.    79 R=ROOTR(I)
  145.       IF(ABS(R)-1.0)81,81,80
  146.    80 R=1.E38
  147.       GO TO 82
  148.    81 J=J+1
  149.    82 CONTINUE
  150.       IF(J-1)83,88,83
  151.    83 IE=2
  152.       GO TO 86
  153. C
  154. C     UNABLE TO COMPUTE R
  155. C
  156.    84 IE=IER
  157.    86 R=1.0E38
  158.       RS=R
  159.       GO TO 100
  160.    88 IF(R-1.0E38)90,83,83
  161. C
  162. C     STANDARD ERROR OF R=0.0
  163. C
  164.    90 RS= SQRT(P1*P2*Q1*Q2)/(Y1*Y2* SQRT(FN))
  165. C
  166.   100 RETURN
  167.       END
  168. 
  169. COF,COF,7,ROOTR,ROOTI,IER)
  170. C
  171.       J=0